Phase diagram of Rydberg atoms with repulsive van der Waals interaction 



O 

(N 

o 

Q 



O. N. Osychenko, 1 G. E. Astrakharchik, 1 Y. Lutsyshyn, 2 El Yu. E. Lozovik, 3 and J. Boronat 1 

1 Departament de Fi'sica i Enginyeria Nuclear, Universitat Politecnica de Catalunya, Campus Nord B4-B5, E-08034 Barcelona, Spain 

2 Institut fiir Physik, Universitat Rostock, D-18051 Rostock, Germany 
3 Institute of Spectroscopy, 142190 Troitsk, Moscow region, Russia 
(Dated: December 30, 2011) 

We report a quantum Monte Carlo calculation of the phase diagram of bosons interacting with a repulsive 
inverse sixth power pair potential, a model for assemblies of Rydberg atoms in the local van der Waals blockade 
regime. The model can be parametrized in terms of just two parameters, the reduced density and temperature. 
Solidification happens to the fee phase. At zero temperature the transition density is found with the diffusion 
Monte Carlo method at density p = 3.9 (h 2 /mCe) 3 ^ 4 , where Ca is the strength of the interaction. The solidi- 
fication curve at non-zero temperature is studied with the path integral Monte Carlo approach and is compared 
with transitions in corresponding harmonic and classical crystals. Relaxation mechanisms are considered in 
relation to present experiments. 

PACS numbers: 32.80.Ee, 64.70.-p, 67.85.-d 
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I. INTRODUCTION 

Rydberg atoms have one electron excited to a high energy 
level. Such atoms exhibit strong and highly tunable interac- 
tions which may have an extraordinarily long range. Optically 
excited from suspended clouds of cold atoms, Rydberg atoms 
interact both between themselves and with the surrounding 
unexcited atoms, resulting in a rich behavior of the Rydberg 
systems. 

Due to the strong interactions, a Rydberg atom shifts the 
levels of nearby atoms sufficiently to prevent their subsequent 
excitation. A large number of studies deal with a local block- 
ade regime. In such a regime a Rydberg atom blocks excita- 
tions in its vicinity, and the atomic clouds may be injected 
with well over 10 3 Rydberg excitations before the existing 
excitations block any further ones iJ-Hl]. Unfortunately, the 
arrangement of the excited atoms in such experiments is not 
directly accessible and has been a subject of intense investi- 
gation. Understanding the ordering of Rydberg atoms may be 
important for interpretation of the experimental results, for ex- 
ample for the antiblockade effect predicted in 14]]. It was also 
suggested that a spatially ordered state may allow for a better 
control over quantum states in such experiments J^. Finally, 
there is an exciting possibility of observing phase transitions 
in these versatile systems, especially to states with long-range 
ordering 

Quantum many-body treatments attempting modeling of 
realistic Rydberg systems have been developed in the past 
HI S IH Ijj Hloll . and were successful in reproducing a num- 
ber of important experimental features Q U0l4l2ll . Due 
to complexity, it is often difficult to consider long-range or- 
der with such calculations. Nonetheless, strong short-range 
spatial correlations between Rydberg atoms were obtained 
in the calculations of Refs. Q Q HHl, as the atoms avoid 
each other due to the blockade. Successful observation of 
the antiblockade effect was also a demonstration of a cre- 



ation of the strong short-range correlations flUl . Possibility 
of long-range ordering (crystallization) of Rydberg atoms was 
recently predicted for systems coupled to specially selected 
chirped laser pulses fl, Il3ll . Ordering was also considered, 
and crystalline phase found, in theoretical calculations of both 
one and two-dimensional optical lattices lll4l4l7ll . Remark- 
able non-commensurate crystalline phases in optical lattices 
emerged in Ref . lfl5tl . 

Given the complex nature of the interactions in the Rydberg 
systems, it is important to know how much of the behavior of 
large assemblies of Rydberg atoms stems directly from the 
pair potential of the interaction between the atoms. For this 
reason we aim to study ordering in the simplest model of the 
Rydberg systems. Because of the large number of Rydberg- 
excited atoms in the experiments, we consider the thermody- 
namic limit. While the results are established in the thermody- 
namic equilibrium, many present experiments with Rydberg 
atoms are too short to reach equilibrium. Thus comparison in 
such cases must be made cautiously. 



n. MODEL AND METHODS 

The dominant interactions in the Rydberg systems are usu- 
ally the Forster-resonant dipole-dipole interactions between 
the excited atoms. It was shown by Walker and Saffman 
lfl8i [l9ll that, given a pair of Rydberg atoms in the same state, 
the interaction will not have zeros as a result of the hyperfine 
structure or alignment of the atoms only if the resonant cou- 
pling is from the s to p states. Furthermore, interactions in the 
s + s — > p + p channels depend only weakly on the hyperfine 
structure of the p states, resulting in a nearly isotropic inter- 
action, to within 1CP 2 . This perhaps in part motivates the use 
of the ns Rydberg states in current experiments Jllr7l fl2ll20ll . 
Neglecting the hyperfine structure, the interaction for this res- 
onance is isotropic and its matrix element is given in terms of 
the Forster defect S as Hill 
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which changes from V = C3/1* 3 to van der Waals' V = 
Cq/t 6 (with Cg = — C3/8) for distances much larger than 
the crossover i?3->6 ~ (— Cq/5) 1 ^ 6 . In the case of a strong 
local blockade, the blockade radius is often larger than the 
crossover distance. In such a case, the excited atoms are more 
likely to be found at distances where the interaction is already 
of the van der Waals type. 

The above arguments motivate the repulsive van der Waals 
model for the Rydberg atoms in the local blockade regime. We 
disregard any energy transfer or interactions with the underly- 
ing gas of the ground-state atoms, and particles are treated 
as spinless bosons in three-dimensional space with the many- 
body Hamiltonian 

Defining the reduced units of length and energy as 

allows us to describe the properties of this model universally 
in terms of just two parameters, the dimensionless density pr$ 
and temperature k B T/E$. The units are selected to satisfy 
Eq = h 2 /mr-Q = Ce/r^. The mass m in Eq. (Q]) is the mass 
of the atom. 

It is important to establish the applicability of the bulk 
phase diagram to finite systems. For a cloud of size R and 
number density p, the tail potential energy per particle can be 
estimated as pCs/R 3 . In order for the phase transition to oc- 
cur at the same parameters in the limited system as in a bulk 
one, it is sufficient that the missing potential energy is much 
smaller than the kinetic energy. In the case Tk B ^> Eq, this 
reduces to R > {/ 'pC 6 /Tk B . When Tk B <C Eq, kinetic 
energy is estimated as H 2 p 2 ^ 3 /m and thus R roipr 3 ) 1 ' 9 . 

III. RESULTS FOR THE PHASE DIAGRAM OF THE 
REPULSIVE VAN DER WAALS GAS 

The phase diagram of the model includes a solid at high 
densities and, at lower densities, a gas phase that Bose- 
condenses at sufficiently low temperatures. To locate these 
phase regions, we employed a number of methods, each suit- 
able in a certain area of the phase diagram. At zero tem- 
perature, the model was treated with diffusion Monte Carlo 
(DMC), a projector method which provides an exact ground- 
state energy for bosonic systems. DMC has been used suc- 
cessfully in the past to calculate the equations of state and lo- 
cate quantum phase transitions for a variety of systems. Tran- 
sitions at non-zero temperature were studied with path inte- 
gral Monte Carlo (PIMC), a first principles method which al- 
lows one to compute the averages of quantum operators by 
summing over the quantum partition function of the system. 
Both DMC and PIMC methods allow one to treat systems with 
several hundred particles under periodic boundary conditions, 



with thermodynamic limits obtained by a suitable extrapola- 
tion. Additionally, classical limits were established with clas- 
sical Monte Carlo calculations. In two regimes the location of 
phase transitions could be expressed in a semi-analytical form. 
In the first case, the transition between superfluid and normal 
gas was expressed in terms of the scattering length of the po- 
tential by means of a known relationship. In the second, the 
solid-to-gas transition was located at low temperatures with 
the harmonic theory. The results are summarized in the phase 
diagram shown in Fig. Q] 

At sufficiently high density, the atoms are expected to form 
a crystalline solid. Summing the potential energy of the per- 
fect lattice structures, we conclude that the preferred symme- 
try is fee. While other structures may be excluded on the en- 
ergetic grounds, the energy of the hep structure is very close 
to that of the fee. The difference between the perfect crystal 
energies, E^ cv — Ef cc = 2 x 1Q~ a {ptq) 2 Eq, is small enough to 
be comparable to or even swamped by the temperature effects 
in present experiments (for example, inRefs. |2lll_2j,|20j]). The 
hep phase is anticipated to be metastable with respect to the 
transition to the fee phase. Zero-point motion and tempera- 
ture effects are expected to keep the fee symmetry preferred 
to hep. In the following, the solid phase is always implied to 
have the fee structure. 

Investigation on the zero-temperature line was done with 
the DMC method j2lll22Tl . For importance sampling in the gas 
phase we used a Jastrow form rii<j{ ex P [ — l/2(&/ r y) 2 ] + 
cxp [— l/2(b/(L — r-y)) 2 ] }, Tij < L/2, for a periodic box of 
size L. The second power in 1/r arises from the cusp condi- 
tion of the scattering problem with the repulsive 1 /r 6 potential 
and is also compatible with the presence of long-wavelength 
phonons l23ll . The parameter b was variationally optimized 
beforehand. The Nosanow-Jastrow wave function was used 
for importance sampling in the solid phase [24,25]. It consists 
of the product of the above Jastrow term and a site-localizing 
Nosanow term Yli cxp [— (r, — U) 2 /2^\, where r.; and ij de- 
note correspondingly the coordinates of the atoms and lattice 
sites, and 7 is the second optimized parameter. The breaking 
of exchange symmetry between particles in the solid affects 
the energy only negligibly l26ll . Within the statistical errors 
of the DMC, results for the energies of the fee and hep lattices 
are indistinguishable and both are lower than the energies de- 
rived using bec configuration. 

While the phase transitions are conventionally reported as a 
function of pressure rather than density, density of the Ryd- 
berg atoms is more accessible and controllable experimen- 
tally. We therefore choose to express the transition locations 
in terms of density, even for the first-order solidification tran- 
sition (in this case one needs to specify the coexistence re- 
gion). We find that the equations of state for the fee solid and 
gas phases cross at the transition density 

p c rl = 3.9 ±0.2, (2) 

expressed in the reduced units with the help of Eq. (Q]). The 
coexistence region of the solid and gas phase at zero tempera- 
ture, determined using the double-tangent Maxwell construc- 
tion, is narrow and is in fact smaller than the above error for 
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the transition density (which arises mostly from the extrapola- 
tion to the thermodynamic limit; calculations were performed 
with up to 256 particles). 

The transition line between solid and gas phases at small 
temperatures can be determined with the harmonic theory 
rf27ll . assuming the Lindemann ratio remains unchanged on 
the transition line. The value of the Lindemann parameter 
at melting may be extracted from the DMC calculations of 
the transition density at zero temperature. The resulting low- 
temperature dependence of the gas-to-solid transition density 
is given by 

yharmonic = C y/ (p - p c ) r* j± , (3) 

where p c is the transition density at zero temperature, Eq. ©, 
and the constant C = 8.0 is determined numerically from the 
dispersion curves of the solid and depends on the interactions 
and geometry of the fee lattice. 

A quantum solid melts at lower temperatures than the clas- 
sical one due to the zero-point motion of the atoms. The 
classical transition was located in the canonical ensemble by 
Metropolis sampling of the Boltzmann factor. As the potential 
energy Cq /r 6 is exactly proportional to the square of the den- 
sity, the transition temperature for the classical system also 
scales exactly as T oc p 2 . We find that 

rclassical = 22 / 3\2_0 (4) 

Kb 

As expected, such scaling removes the Planck constant from 
the classical transition temperature, which in fact simplifies to 

T cla S sical fcB = .22p 2 C* 6 . 

To fully account for quantum effects, the gas-to-solid tran- 
sition at T ^ was also located with PIMC calculations. 
We used decomposition of the action operator that is accu- 
rate beyond the fourth order 12811 . For details of the method 
and implementation, see Ref. 12911 . The transition was located 
by observing melting or solidification while working in the 
canonical ensemble, beginning with configurations of atoms 
placed on a randomly distorted lattice. Used in this way, the 
calculations determine a range in which the transition density 
is located. PIMC results confirm the validity of the harmonic 
approximation at low temperatures. At higher temperatures 
the transition density follows the classical melting curve (0). 

The above results establish the solidification transition of 
the repulsive van der Waals model. Additionally, the dynamic 
nature of the Rydberg gas raises a possibility for the spatial 
ordering to be induced kinetically, as the combination of de- 
cay and strong blockade will favor supplanting excitations to 
be equidistant from their immediate neighbors. We modelled 
such a process and observed that replacement of decaying ex- 
citations in the local blockade regime indeed creates a short- 
distance order, but not a true long-distance crystalline order- 
ing. These finding are consistent with much more elaborate 
dynamic models of Ref s. HUE!. 

At low temperature, the gas phase of the model is expected 
to form a Bose-Einstein condensate (BEC). Transition be- 
tween the BEC and normal gas phases at low densities lies 
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FIG. 1. (Color online) Phase diagram of the repulsive C(,/r b inter- 
action, scaled to units given in Eq. iQJ. Location of the gas-to-solid 
transition at zero temperature, determined with DMC, is shown with 
an open blue square on the T = axis. Dashed blue line shows gas- 
to-solid transition as found with harmonic theory [Eq. l[3}]. Solid 
blue line shows classical gas-to-solid transition [Eq. ©]■ Solid blue 
squares show location of the gas-to-solid transition as determined 
with PIMC. Red short-dashed line marks the Bose-Einstein con- 
densation of the ideal gas [Eq. {5J]. Filled red bullets show Bose- 
Einstein condensation temperatures found with PIMC. 



slightly above the ideal Bose gas condensation temperature, 

(3 \ 2/3 j-> 

due to the repulsive interaction between particles & The 
correction is governed by the scattering length of the po- 
tential a s , which can be found to be equal to a s = 
2r(3/4)/r(l/4)r = 0.676. ..r . The transition temper- 
ature is then given by Tbec = (l + casp 1 ^ 3 ), where c 
is a positive constant of the order of unity (for details, see 
Ref. 1I30I1 and references therein). In the present case this ex- 
pression is only valid at very low densities (one needs to sat- 
isfy at least pr^ < 5 x 10~ 2 to make the description in terms 
of the zero-momentum scattering length meaningful), where 
the magnitude of the correction is not significant. 

At higher densities the BEC-to-normal gas transition is no 
longer universal and depends on the form of the potential. 
We determined the location of this second-order transition 
with the PIMC method by calculating the superfluid transition 
from the winding number estimator J3lll . The PIMC calcu- 
lations show that, at higher densities, the interactions deplete 
the condensate and the transition temperature is lower than 
for the ideal Bose gas. Combining the PIMC results, the re- 
gion in which the triple point is located was determined as 
4.5 < T/(E a /k B ) < 6.5 and 4 < pr^ < 5, which we con- 
sider sufficiently narrow for practical considerations. 
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IV. COMPARISON WITH EXPERIMENTAL CONDITIONS 

Because the interaction constant C% enters the reduced units 
[Eq. (Q])], the effective temperature and density can be var- 
ied over many orders of magnitude. Most of the present ex- 
periments are deeply in the "classical" region of the phase 
diagram (Fig. [T]). As an example, we consider the condi- 
tions of the experiments presented in Ref. lEoll . For the ex- 
citation with 170 ns laser pulses, the system parameters at 
4 pK are T/(E Q /k B ) w 33 x 10 5 and prl w 1.9 x 10 3 , 
which in fact correspond to the gas phase of the equilibrium 
phase diagram. For 320 ns excitation pulses and T = 1 pK, 
T/(E /k B ) « 8.2 x 10 5 and prl w 7.4 x 10 3 , well below 
the gas-to-solid transition. Therefore, the achievable temper- 
ature and density are already in the range suitable for investi- 
gating the equilibrium phase diagram. Increasing the excita- 
tion number increases the interaction constant Cq and moves 
the system deeper into the classical regime where the gas and 
solid phases are separated by the simple condition of Eq. (@). 
The quantum regime of the phase diagram may be accessed 
by decreasing the excitation numbers or increasing the Ftirster 
defect 5. 

Whether Rydberg atoms in actual experiments will reach 
or even approach an equilibrium phase depends on their life- 
time, the experiment duration and availability of relaxation 
mechanisms. Because of the short lifetimes of the Rydberg 
states, most current experiments are performed on such short 
timescales as to make the thermal motion negligible. It is 
therefore said that the experiments are performed with Ryd- 
berg excitations of a frozen gas. If the experiments are ex- 
tended closer to the currently achievable lifetimes of the Ryd- 
berg states, which can be as large as 100 ps HH], some 
degree of thermal equilibration will already be achieved. Be- 
sides the thermal motion there are, however, at least two other 
kinds of motion that may need to be considered. The first one 
is the motion of the excited atoms due to the strong forces be- 
tween them. The characteristic timescale associated with such 
a motion is the time that it takes for a Rydberg atom to travel 
the mean distance between Rydberg atoms. Given the mean 
distance £ « p -1 / 3 and the imbalance force of the order of 
Cg /£ 6 , this time is given by 



^ballistic 



C 6P ' 



8/3 ' 



(6) 



For a small fixed number of Rydberg excitations time © de- 
creases rapidly with the excitation number n as n~ n / 2 ; as 

— 1/2 

the local blockade is reached, p oc C 6 , and tbaiiistic in- 
stead grows as n 11 / 6 . For example, Rydberg systems created 
by the 1.970 ps pulses from lpK gas in the experiment of 
Ref. rf20ll have tbaiiistic ~ 12 ps. For the setup of Ref. 01, 



^ballistic ~ 60 ps while the clouds could be successfully stud- 
ied for as long as 20 ps. Collisional ionization and heating 
could potentially hamper such relaxation 13311 . 

V. DISCUSSION AND CONCLUSIONS 

As an example, consider the possibility of exploring the su- 
persolidity in Rydberg systems. Ground state atoms dressed 
in Rydberg states exhibit weak van der Waals interactions at 
large distances, as described in H34II. S upersoliditv of such 
atoms was already investigated 11341. 13511 . Consider instead a 
scenario in which the gas is additionally allowed to have a lat- 
tice of Rydberg excitations. Such a lattice would in turn im- 
pose weak but long-range spatial correlations onto the ground- 
state atoms. At the same time, the ground-state atoms may 
be placed in the BEC regime 12011 . However, it may be im- 
possible to identify which of the atoms was excited within a 
certain proximity, as was demonstrated, for example, by the 
superatom analysis of the experimental results in Refs. 01 El- 
(However, motion will lead to dephasing of this state 0611 .) If 
the atoms are indeed prepared in such a mixed state, com- 
bining the ground \g) and excited |e) states as \g) + a\e), 
Ng 1 <C a. <C 1, then both the lattice-forming and the 
BEC components are indistinguishable and may be said to be 
formed by the same atoms. Therefore, such a system would 
consist of particles which would simultaneously break trans- 
lational symmetry and possess off-diagonal long-range order, 
an epitomic realization of supersolid. While our model does 
not include the light field, the conditions for the phases of 
both excited- and ground-state atoms may be immediately ex- 
tracted from Fig. Q] just with different reduced units for the 
two species. 

In conclusion, it is possible to parametrize a model with 
isotropic van der Waals interactions into a universal phase 
diagram. We have characterized the phase diagram of Ryd- 
berg atoms by considering a model of bosons with repulsive 
van der Waals interaction, and determined solidification and 
Bose-Einstein condensation conditions. Relaxation mecha- 
nisms other than thermal motion should be considered if one 
considers Rydberg systems on time scales of several tenth of 
microseconds. Finally, it is worth mentioning that interactions 
between Rydberg excitations open a possibility of new super- 
solid scenarios. 
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